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ABSTRACT 

We obtain self-similar solutions that describe the gravitational collapse of 
nonrotating, isothermal, magnetic molecular cloud cores. We use simplifying 
assumptions but explicitly include the induction equation, and the semianalytic 
solutions we derive are the first to account for the effects of ambipolar diffusion 
following the formation of a central point mass. Our results demonstrate that, 
after the protostar first forms, ambipolar diffusion causes the magnetic flux 
to decouple in a growing region around the center. The decoupled field lines 
remain approximately stationary and drive a hydromagnetic C-shock that moves 
outward at a fraction of the speed of sound (typically a few tenths of a kilometer 
per second), reaching a distance of a few thousand AU at the end of the main 
accretion phase for a solar- mass star. We also show that, in the absence of field 
diffusivity, a contracting core will not give rise to a shock if, as is likely to be 
the case, the inflow speed near the origin is nonzero at the time of point-mass 
formation. Although the evolution of realistic molecular cloud cores will not be 
exactly self similar, our results reproduce the main qualitative features found in 
detailed core-collapse simulations (Ciolek & Konigl 1998). 

Subject headings: accretion, accretion disks — diffusion — ISMxlouds 
— ISM: magnetic fields — MHD — stars: formation 
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1. Introduction 

Low-mass stars are generally believed to form as a result of the gravitational collapse 
of molecular cloud cores. The cores are initially supported by thermal and magnetic 
forces, but because of ambipolar diffusion (the drift of ions, to which the magnetic field 
lines are attached, relative to the dominant neutral gas component), they gradually lose 
their magnetic support and eventually collapse after becoming "supercritical" (see, e.g., 
Mouschovias 1987 for a review) .0 The most detailed numerical treatments to date of the 
problem of the ambipolar diffusion-initiated formation of supercritical cores and the early 
stages (prior to point mass formation) of their subsequent dynamical collapse have been 
presented by Mouschovias and collaborators (Fiedler & Mouschovias 1992, 1993; Ciolek & 
Mouschovias 1993, 1994, 1995, hereafter CM93, CM94, CM95; Basu & Mouschovias 1994, 
1995a, 1995b, hereafter BM94, BM95a,b). Because the timescale for core formation is much 
longer than the timescale for dynamical collapse, special numerical techniques had to be 
employed in these calculations. The simulations were terminated when the central densities 
reached ~ 10^° cm~^ and the underlying assumptions of isothermality (e.g., Gaustad 1963) 
and flux freezing onto the ions (e.g., Pneuman & Mitchell 1965) broke down. These 
calculations were nevertheless able to demonstrate that supercritical cores begin to collapse 
dynamically before a point mass (i.e., a protostar) appears at the origin. 

The dynamical evolution of supercritical cores after their formation has been studied 
by many researchers. Solutions exist for the collapse of nonrotating, self-gravitating spheres 
without thermal support (Henriksen 1994), self-gravitating spheres with thermal support 
(Penston 1969; Larson 1969; Shu 1977; Hunter 1977; Boss & Black 1982; Whitworth & 
Summers 1985; Foster & Chevalier 1993) as well as with a combined thermal and isotropic 
magnetic pressure support (Chiueh & Chou 1994), and self-gravitating disks with thermal 
support (Narita, Hayashi, & Miyama 1984; Matsumoto, Hanawa, & Nakamura 1997) and 
also with ordered, frozen-in magnetic fields (Nakamura, Hanawa & Nakano 1995; Li & 
Shu 1997, hereafter LS). In order to choose a particular solution for a given problem, one 
needs to know the properties of the supercritical core at the time of its formation. This 
information, however, can only be gleaned from a study of the preceding, quasi-static 
evolution of the core under the influence of ambipolar diffusion. Although different 
assumptions about the initial state of the core yield solutions that are qualitatively similar 
in their gross behavior (the core collapses with near free-fall speeds and a point mass 
eventualy forms at the center), the solutions do differ in such important details as the 



^In this paper we reserve the term "core" for the high-density central region of a molecular cloud and do 
not apply it to the point mass that forms from the collapse of such a core. 
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accretion rate onto the central point mass and the formation (or absence) of shocks. 

The well-known examples of the Larson- Penston (1969) and Shu (1977) similarity 
solutions in fact represent two extremes of a whole continuum of self-similar collapse 
solutions specified by a cloud's initial configuration and the conditions at its boundary 
(Hunter 1977; Whitworth & Summers 1985; see also Chiueh & Chou 1994 for a 

generalization to the case of an isotropic internal magnetic pressure). The Larson- Penston 
(1969) solution is characterized by a spatially uniform, supersonic (at ~ 3.3 times the 
isothermal speed of sound C) infall speed and an inverse-square dependence of the density 
p on the radius r at the instant of point-mass formation (PMF); the mass accretion rate 
at the center is ~ 29 C^/G (where G is the gravitational constant) at that instant and 
increases to ~ 47 C^/G immediately after PMF. Numerical simulations of the collapse of 
nonmagnetic isothermal spheres (Hunter 1977; Foster & Chevalier 1993) have indicated 
that this solution provides a good approximation to the conditions near the center at 
the PMF epoch for clouds that are initially near a marginally stable equilibrium. The 
Shu (1977) solution strictly applies only to the post-PMF evolutionary phase: it consists 
of an inner free-fall region and a hydrostatic outer envelope that are separated by an 
outward-propagating (at a speed C relative to the gas) expansion wave. The envelope 
corresponds to a singular isothermal sphere (p oc r~^) and the mass accretion rate onto the 
center is ~ 1 C^/G. In applying this solution to real systems, it was proposed to identify 
the initial core configuration at the end of the quasi-static ambipolar-diffusion phase with a 
singular isothermal sphere (or, more generally, a toroid) at the instant of PMF (e.g., Shu, 
Adams, & Lizano 1987; Li & Shu 1996). However, as we noted above, the conclusion from 
detailed numerical simulations has been that the dynamical phase of core collapse generally 
commences well before the PMF epoch, so that the innermost region is not well represented 
by a quasi-static solution at the time of point-mass formation. 

Another interesting effect that depends on the specific choice of initial conditions and 

on the detailed physical properties of the collapsing core is the formation (or absence) 
of shocks (e.g., Tsai & Hsu 1995; LS). For example, LS discovered that when, instead 
of a spherical core, one considers the collapse of a flattened disk, the expansion wave of 
Shu (1977) becomes a shock. As we show in this paper, when one takes proper account of 
the fact that supercritical cores collapse dynamically before a point mass first forms at the 
origin, that shock disappears. Nevertheless, a physical basis for the formation of shocks in 
collapsing magnetized molecular cloud cores has been discussed by Li & McKee (1996), who 
argued that a hydromagnetic C-shock will appear as a result of the outward diffusion of 
inwardly advected magnetic flux. The existence of such a shock has been confirmed in the 
numerical simulations of Ciolek & Konigl (1998, hereafter CK), and it is, in fact, a salient 
feature of the semianalytic solutions derived in this paper. 
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The aim of the present work is to clarify the effects of ambipolar diffusion in 
dynamically collapsing supercritical cores. Toward this goal, we construct semianalytic, 
time-dependent similarity solutions of gravitationally contracting, magnetized, isothermal 
disks. Although the evolution of real molecular cloud cores is not expected to be exactly self 
similar, we demonstrate, through a comparison with the detailed numerical simulations of 
CK, that our solutions capture the main traits exhibited by the latter calculations. Based 
on an analogous comparison with the results of numerical simulations, Basu (1997) showed 
that a self-similar scaling describes the pre-PMF evolution in the innermost flux tubes 
of collapsing supercritical cores quite well. To complement his study, we concentrate in 
this paper on the post-PMF evolutionary phase. Our approach differs, however, from that 
of Basu (1997) in that we explicitly solve the induction equation, whereas he accounted 
for the effects of ambipolar diffusion only in a phenomenological manner.^ In fact, the 
solutions that we derive, while involving various simplifications, are nevertheless the first to 
consistently incorporate ambipolar diffusion into a self-similar representation of the collapse 
of a magnetized cloud core. We formulate the problem in §2, present our solutions in §3, 
and discuss the results in §4. Our conclusions are summarized in §5. 

2. Mathematical Formulation and Approximations 

The problem of the self-initiated formation and contraction of cloud cores due to 
ambipolar diffusion in axisymmetric, self-gravitating, isothermal, magnetic molecular clouds 
was formulated in detail in CM93. In the presence of an ordered, large-scale magnetic field 
the contracting cloud core assumes a disk-like configuration on account of the fact that the 
magnetic stresses inhibit motion normal to the field lines (see Fiedler & Mouschovias 1993). 
Along the field lines it is a good approximation to assume that the cloud is at all times in 
hydrostatic equilibrium, with thermal pressure providing support against vertical gravity. 
For simplicity, we neglect both the pressure support provided by internal hydromagnetic 
waves (i.e., "turbulence") and the vertical squeezing of the core by the external thermal 



■^Our work is thus also distinguished from that of Safier, McKee, & Stahler (1997), who studied the effects 
of ambipolar diffusion in the spherically symmetric, quasi-static limit without explicitly solving the induction 
equation. 

"^The effect of weak magnetic fields on a dynamically collapsing core in the presence of ambipolar diffusion 
was previously investigated by Galli & Shu (1993a), who carried out a perturbation expansion of the 
(nonmagnetic) spherical similarity solution of Shu (1977). As was already noted and discussed by Li & 
McKee (1996), the semianalytic solution derived in that paper, as well as the associated numerical calculation 
in Galli h Shu (1993b), did not uncover the existence of a flux diffusion-driven shock. 
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and magnetic pressure (including, in particular, the squeezing induced by the radial field 
component at the disk surface; see eq. [9] in CK). For a self-gravitating disk with possibly 
a point mass at the center, the hydrostatic balance equation can thus be written as 

pC' = ^G.V^^. (1) 

where r is the cylindrical radius, a is the surface density, p is the gas density, h = a/2p is 
the vertical scale height (assumed to be <C r), C = (ksT /rrinY^'^ is the isothermal speed of 
sound (with ks being the Boltzmann constant, T the temperature, and m„ the mean mass 
of a gas particle), and Mc is the mass of the star that forms at the center. Equation (|l]) 
implies that 

_ / rC{GMj2r)-^/^ near the central point mass , , 

I C'^lTcGa]^^ far from the central point mass . 

The collapse of the core can be described using the mass and (radial) momentum 
conservation equations, which, in cylindrical coordinates (r, 0, z) , read 

da 1 d , , 

a^ + ;g;(''-) = o (3) 

and 

dt dr ^ dr 2'n-a \ ^''^ dr J ' 

respectively. Here i?r,s is the radial component of the magnetic field at the surface of the 
thin disk (by symmetry, Bj. = ai the midplane z = 0, and dBr/dz ^ Br^s/h), u is the 
radial velocity component, and is the radial component of the gravitational acceleration. 
As is well known, the calculation of Qr for a self-gravitating thin disk is in general nontrivial. 
The above equations are supplemented by the induction equation for the magnetic field (see 
derivation below). 



1/2 / 

UB,. + I '-^^ 1 Br,s 





(5) 



where 

r/ ^ r„,/(87rGp)-^/2 ^ 0.7 Clif , (6) 

the ratio of the mean collision time r„j of a neutral particle with a sea of ions and the 
approximate free-fall time due to the disk self-gravity, is a measure of the efficiency of 
ambipolar diffusion in the disk, and 

(T) 
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is the dimensionless local mass-to-flux ratio. In equation \l/(r, t) = 2TTr'Bz{r' ,t)dr' 
is the magnetic flux threading the disk interior to the radius r at time t, whereas in 
equation (|^) ^ = 10^^^,^_i7 is the cosmic-ray ionization rate per hydrogen nucleus (e.g., 
Watson 1976). In deriving equation we have assumed that the disk is weakly ionized 
and that the magnetic fleld is "frozen" into the ions (valid for disk densities <^ 10^^ cm~^). 
The flrst assumption implies 

Pi < Pn , (8) 

where the subscripts i and n refer, respectively, to the ions and neutrals. Under this 
assumption, one can identify p in the preceding equations with p„, and u with m„. The 
speeds Un and Ui are not, in general, equal, and their difference ud = iui — Un) is referred to 
as the (radial) ion-neutral drift speed. The second assumption can be expressed by writing 
the flux conservation equation in the form 

— = -2nrUiB^ = -2nr{un + Ud)B^ . (9) 

Magnetic forces that act on the ions are transferred to the neutrals through a collisional 
drag force, 

pv^ _ fOB^ dBr\ , . 

Tm ~ Ait[ dr dz J - ^ ' 

For a disk that consists of a molecular hydrogen gas with a helium number density that 
is 10% of the density of hydrogen nuclei, the neutral-ion mean collision time is given by 
Tni = 1.4(mi + mH2)/[Pi < (^ui >iH2]y where < aw >iH2= 1-69 x 10~^ cm^ s~^ is the average 
collisional rate between ions and hydrogen molecules (e.g., McDaniel & Mason 1973). By 
balancing ionization by cosmic rays with dissociative recombination, one can express the 
ion number density as rii ^ (^n/fldr)^/^ (Elme green 1979), where ajr ~ 10 ^ cm'^ s ^ is the 
electron dissociative recombination rate (e.g., Dalgarno 1987) and n is the neutral particle 
number density. Although this expression (which was used in deriving the dependence of rj 
on ^ in eq. 0) is strictly valid only for densities n ^ 10^ cm~^ (at densities S> 10^ cm~^, 
rii ^ const] see Figs. 2c and 4c in CM94, and Figs. 1 - 3 in Ciolek & Mouschovias 1998), 
we nevertheless adopt it in this paper for all disk densities since it allows us to obtain a 
self-similar set of equations that can be directly integrated. Equation (|^) follows from 
combining equations (|^) and (p^ and taking the ion mass to be nii = 29 uih- Equation (|I0|) 
brings out the well-known fact that, in order to support a nonuniform magnetic fleld 
distribution in the disk, the ion-neutral drift speed has to be nonzero. We return to this 
basic point in §3.1. 

Given the distributions a{r, 0), u{r, 0), and Bz{r, 0) at some initial time t = 0, one can 
in principle integrate equations (0)-(||) provided that gr{r,t) and Br^s{r,t) can be expressed 
as functions of a and B^, respectively. In the limit of a thin disk with negligible mass 
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outside, one can obtain the radial gravitational acceleration at radius r by adding the 
contributions from rings of disk material at all radii r', 



CO 



G 

g^{r,t) = — - I 2nr'a{r',t)n{r'/r)dr' . (11) 



In this equation, 

1?(Y\ ^ ['^ (1-Xcos0)d0 d/C 

^^^^ = 2^70 (l + X^-2Xcos0)3/^=^^'')+''dX' ^''^ 

where 

m) = ^^(TTx) ^^^^^'^ ^ ^^'^ ' ^^^^ 

with i^(x) being the complete elliptic integral of the first kind (see CM93). Figure 1 
shows a plot of 'R-{X) obtained using approximate formulas for K{x) (e.g., Abramowitz & 
Stegun 1965, p. 591). It is seen that contributions to the gravitational acceleration at a 
radius r arise from disk material both interior and exterior to r. There exist, however, two 
limits in which the simple expression 

holds true, where M{r,t) = Jq 27Tr'a{r' ,t)dr' is the mass interior to the radius r at time 
t: (a) when a oc r~^, and (b) when the gravitational pull of a point mass at the center 
dominates over the disk self-gravity. As we demonstrate below, equation (|1^ is a very good 
approximation to the true value of Qr during the post-PMF phase of the disk evolution. For 
future reference, we note here that the mass conservation equation (|^) can be rewritten in 
terms of M(r, t) in the form 

dM DM ^ 

What about -B^.s? Assuming that the medium surrounding the disk is current-free, 
the magnetic field can be expressed as the gradient of a potential: B = V$. Now, since 
V ■ B = 0, this scalar potential satisfies Laplace's equation 

= (16) 

with Neumann boundary conditions 

9$ J Bz at the surface of the disk {z ~ 0) 
dz 1-82,00 at infinity , 

where -82,00 is the uniform and constant large-scale magnetic field that is frequently observed 
to thread molecular clouds on scales that are large in comparison with their (flattened) inner 
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cores (e.g., Hildebrand, Dragovan, & Novak 1984; Novak, Predmore, & Goldsmith 1990; 
Kane et al. 1993). Therefore, in the ideal case of a current-free medium outside a thin disk, 
the function $ — zB^^oo satisfies the same equation and (Neumann) boundary conditions as 
the gravitational potential V (where g = —W, \/^V = 0, dV/dz = 27cGa at z = 0, and 
dV/dz = at infinity). Hence, by direct analogy with the expression for Qr in equation ([TT|) , 
one can write 



2=0 



1 

- r\B,{r',t)-B,^^)n{r'/r)dr' 
Jo 

1 

— / r'BJr',t)n(r'/r)dr' , 
Jo 



(17) 



'1^ 



where the last approximation is valid when B^ at the disk surface is ^ -82,00 (see CM93, 
CM94, CM95, BM94, and BM95a,b). When one considers only the inner parts of the 
collapsing disk, as in the present analysis, one can neglect -62,00 -Q Pursuing the analogy 
with gravity even further, one sees that it is possible to use the simple expression 



Br 



27rr2 



(19) 



when either (a) B^ oc r^^ or (b) the magnetic field advected to the center can be represented 
by a split monopole that dominates over the disk magnetic field. We will demonstrate that 
either (a) or (b) are applicable in a collapsing disk following the formation of a central point 
mass, so that equation (|T9|) provides a good approximation to B^^^ during that phase. 



Self-Similar Solutions 



3.1. The "Pivotal" State 



Although one can in principle integrate equations (^-(0) for any set of physically 
reasonable initial conditions [expressed by the functions cr(r, 0), u{r,0), and Bz{r,0)], one 
particular set of initial distributions makes it possible to obtain self-similar solutions that 



'"'When the strong inequality 3> -62,00 no longer holds, as is the case in the outer parts of the 
supercritical core, the expressions in equations ( pl| ) and ( p^ ) differ, and magnetic tension can indeed 
"overwhelm" self-gravity. This should clarify the issue raised in the footnote on p. 248 of LS. 
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compare very favorably with detailed numerical simulations.^ Specifically, the only initial 

surface-density and magnetic-field distributions that are consistent with a self-similar 

evolution are ^ 

a(r,0) (X BJr,0) = Brsir,0) oc - . (20) 

r 

For these distributions the initial mass-to-flux ratio is constant throughout the disk:[] 

a(r,0) _ Xo 



const. (21) 

As it turns out, the scalings given by equation (|20| ) are indeed representative of the central 
region of a collapsing core at the time of point-mass formation (see CM94, CM95, BM94, 
BM95a,b, and Basu 1997). We therefore associate that instant with the time t = 0, and, 
adopting the nomenclature introduced by Li & Shu (1996), we denote the corresponding 
disk configuration as the "pivotal" state of the model. 

We emphasize once again the point first made in §1 that the pivotal state is not the 
end state of the sequence of subcritical quasi-static equilibria that slowly evolve due to 
ambipolar diffusion. | In fact, our pivotal state represents a dynamically collapsing disk 
with 

M(r, 0) ^ , u{r, 0) ^ . (22) 

A simple physical argument clarifies why the self-similar pivotal state cannot correspond 
to hydrostatic equilibrium when ambipolar diffusion is self-consistently taken into account. 
Self-similarity requires a nonuniform distribution of -B^ at t = (eq. [Q), which, in order 
to be supported, requires a nonzero (and positive) drift speed ur, between ions and neutrals 
(eq. [l^l); otherwise Br^s{r, 0) would have been zero and Bz{r, 0) would have been uniform, 
which is clearly not the case. During the initial phase of core collapse (t < 0) the magnetic 
flux either remains fixed in space or else is advected inward, i.e. Ui{r,t < 0) ^ 0. Hence 

\u{r,0)\ = \u^{r,0)-UD\;>UD = C^(l + ^j , (23) 



^Note that, although the numerical models exhibit a self-similar behavior at late times irrespective of 
the detailed early-time {t — s- —oo) conditions, only one particular set of distributions corresponds to exact 
self-similarity at t = 0. 

^LS refer to disks in which the mass-to-flux ratio is a spatial constant as being isopedic. 

^It is interesting to note, however, that our t — surface-density and magnetic field distributions (eq. 
pO{ ) are the same as those of the static pivotal state in the LS similarity solution, which corresponds to a 
singular isothermal disk (see Li & Shu 1996). 
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where the expression for ud follows from equations (Q), (p!0[), (pO|), and (0). One can obtain 
an order-of-magnitude estimate for M(r, 0) by using also equation (|^) to approximate 
(j(r, 0) by its value in equilibrium, 

, , C2 (\l + 2\ 3r]C , , , 3r) 

°' 2S> (AbT j • (A5T2) • ^''<'-- «' G (AFT) ' 

Taking as representative values Ao = 2.9 and ^_i7 = 5 (see §3.4), one infers 

/ J' \ 3/2 

|M(r, 0) I ;^ O.IC and M(r, 0) ^ 0.2 — — M© per 10^ yr , (25) 

\10 K/ 

where T is the temperature. The numerical simulations of CK suggest that, in fact, 
\u{r, 0)1 ^ ud, and that M(r, 0) is of the order of a few solar masses per 10^ years.0 



3.2. Dimensionless Equations 



m(x) = / a{x')x'dx' , m(x) = —axv , '?/'(x) = / 62(x')x'dx' 



We start by introducing a similarity variable x and dependent nondimensional variables 

h{x), a{x), v{x), g{x), m{x), m{x), b(x), and ^{x): 

X = r/Ct , (26) 

h = Ct h{x) , (T(r, t) = [C/{2iTGt)] a{x) , u{r, t) = C v{x) , = (C/t) c/(a;) , (27) 

M(r, t) = (C^t/G) m(x) , M(r, t) = (CVG) m(x) , (28) 

B = {C/[G^'h]) b(x) , ^(r, t) = i){x) , (29) 

where 

i"r. t"T. 

(30) 

(31) 
(32) 

(33) 

(34) 



In these variables, equations (H)- 

da a 
dx 1 — fx — t;) 



and (|l|) can be rewritten as 

(x — vY 



X 



dv . d In a x — v 

— = (x-v)— — + 

dx dx X 



h 



x{x - v)b;, - XT] \— j I — I 16^ 



dx , 



ax 



2mr 



-1+1 + 



orric 
a^x^ 



1/2- 



^Note that, for a 10 K i/2 gas with an interstellar He abundance, irin = 2.33 a.m.u. and C^/G — 
1.58 M0 per 10^ yr. 



- 11 - 



where 

g = / x'a(x')}C(x' /x)dx' , K s = ^ / x'bJx')IC(x'/x)dx' , (35) 

x^ Jo ' x^ Jo 

and rric = m(0). These equations are integrated subject to the initial conditions 

a— > — , bz ^ — , V ^ Vo , m ^ Avo = rho as x ^ +00 . (36) 
x Xo 

Since the value (mc) of the central point mass is not known a priori, the solution requires a 
numerical iteration to calculate h. 



The complicated expressions for g and br,s in equation (|35|) can in principle be solved 
by means of an iterative procedure. As noted by LS, it is natural to start such an iteration 
using the monopole terms 

9---, Ks ~ ^ , 37 
x^ x^ 

where, by virtue of equation (|15|), m(x) satisfies 

m = x{x — v)a . (38) 

With these substitutions one can solve for a{x) and 62 (x), and, in turn, use the latter 
distributions to improve the estimates for g{x) and br^s{x)- This procedure can be repeated 
until a good convergence is attained. As it turns out, the behavior of the disk at t > is 
already very well described by the monopole approximation, so, in what follows, we only 
use the monopole terms in deriving our solutions. 



The system of equations (pl|)- (|3^) contains a singular line at the locus of points in the 
(x, —v) plane where 

{x-vf = 1 , (39) 

which corresponds to the sonic line. This result should be contrasted with the criticality 
condition for ideal MHD, which involves the magnetosonic speed instead of just the thermal 
sound speed (see §3.3). The mathematical reason for this difference is that magnetic 
diffusivity introduces derivatives of higher order than in ideal MHD (e.g., Ferreira & 
Pelletier 1995). There is also an interesting physical explanation of this result. In a 
magnetohydrodynamic flow, information propagates via magnetosonic and sonic waves. 
The critical surfaces that appear in steady-state problems can be thought of as being the 
relics of the time-dependent problem in which the various waves had sufficient time to 
propagate to the exterior boundaries and communicate the information associated with 
those boundaries to the whole flow (Blandford & Payne 1982). When diffusivity is present, 
the only waves that can survive propagation to and from infinity are the sonic waves (since 
the magnetosonic waves simply dissipate away): this is the physical reason why only the 
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sonic critical point appears in steady-state, diffusive MHD. The collapse problem that we 
have formulated, although time dependent, resembles a steady-state problem in this respect 
because of the restrictive assumption of self-similarity, which effectively combines time and 
space into a single variable. It therefore also involves the establishment of a critical surface, 
whose nature is determined by the above argument. 



3.3. The Flux-Frozen Case 

We first consider the ideal-MHD case f] = 0, with Vq ^ 0. In this limit, one can rewrite 
equation (|55|) (using eqs. and as 



din _ ding ^^^^ 
da; dx 

which introduces a singular line in the {x, —v) plane given by 

{x-vf = 1 + 2X-^ (41) 

(where, under the assumed ideal-MHD conditions, the mass-to- flux ratio A = a/62 is a 
constant). This singular line is different from the one given by equation ( |5^ for the 
nonideal case. In physical units, equation (|4l|) can be rewritten as r/t — u = ±(C^ -|- u\Y^'^, 
where ua = -Bz/(47r[7rG'cr^/2C^])^/^ and (C^ + f^)^^^ are, respectively, the effective Alfven 
and magnetosonic speeds in the disk.|^ The equations now describe the t > collapse of an 
isothermal and isopedic disk with X{x) = Xo =const. The solution for the representative 
parameter set Xo = 2.9, —Vo = 1, and rho = 3 (see §3.4) is shown in Figures 2a-2d. It is 
seen that, as the collapse progresses, mass accumulates at the origin at a rate 

/ 7^ \ 3/2 

= 9.6 — - Mq per 10^ yr , (42) 



,10 K. 

corresponding to rric = 6.1 in the expression 

M,(r = 0+, t) = (CH/G) rUc (43) 

for the central point mass. Note that, as discussed in connection with previous collapse 
calculations (see § 1), there is a significant and rapid increase (rhc = 2.0 rho) in the accretion 



^°When the density is given by p = 7rG(T^/2C^ (see eq. ua is equal to the usual Alfven speed in the 
disk. The expression (C^ + given by Shu & Li (1997, footnote 2) for the effective magnetosonic 

speed is evidently an error. 
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rate following point-mass formation. Near the origin, the flow proceeds in "diluted" free 
fall, 

-v= [2m,(l - A;2)/x]i/2 , a = K/2(l - \-/)xY''' , m = . (44) 
The magnetic field advected to the center assumes a split-monopole topology, 

0^ = — (XX ' , hr^s = T 2 ^ ■ 

These results, obtained using the monopole approximation (eq. |3^), change little when 
one corrects for the deviations from the exact expressions for g and 6^,5- They differ from 
the ideal-MHD results presented in LS in that the collapse occurs without discontinuities in 
the flow parameters or their derivatives. The robustness of this conclusion is manifested 
by the fact that the solution curve in the (x, —v) plane never approaches too closely to 
the singular line. As we noted in §1, the initial state adopted by LS is characterized by 
fo = 0: as a result, their integration from x = +oo toward x = encounters the singular 
line. In order to form a point mass at the origin, they then need to connect to one of 
the "minus" solutions defined in Shu (1977). When they perform the integration in the 
monopole approximation, the connection takes place with continuous flow parameters (but 
discontinuous first derivatives) at the point where the singular line crosses the axis —v = 0. 
However, when they use more accurate expressions for Qr and Br^s, the curve intersects the 
singular line a little bit below the axis —v = 0, where (according to Fig. 2 in Shu 1977) no 
"minus" solutions exist. In this case the connection to a "minus" solution can be achieved 
only through a discontinuous jump across the singular line, and this is the weak shock that 
LS invoke. Our solution demonstrates that such a shock is not required when Vo 7^ 0. This 
result is confirmed by the fiux-frozen collapse model presented in §3.3 of CK. (As discussed 
in § 1, numerical simulations of nonmagnetic collapse have shown that the Larson-Penston 
1969 solution, which also has Vo 0, provides an accurate approximation to the physical 
state at the instant of PMF in unmagnetized clouds.) 



3.4. The Ambipolar-DifFusion Shock 

We now proceed to derive the solution for the realistic, nonideal case (i.e., f] ^ for 
X > 0). We continue to employ the monopole approximation for the gravity and magnetic 
field, which we subsequently justify both by verifying that the correction terms remain 
small and by comparing the results with the CK numerical simulations. One immediately 
realizes that the split-monopole field topology (eq. |^) that characterizes the ideal-MHD 
solution near the origin cannot apply in the presence of ambipolar diffusion. This is because, 
for the split monopole, br^s oc ^ bz (x x"^/^ as a; ^ O'^, so the diffusive term in the 
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induction equation (eq. 0) greatly exceeds the advective term near the origin. In practice, 
the outward diffusion prevents the formation of the spht monopole field configuration in 
the first place, and a different solution is established near the origin. This solution is 
characterized by the magnetic flux being left behind while matter continues to fall into the 
center and forms there a point mass that grows with time according to equation (^3|). The 
flow consequently turns into an "undiluted" free fall at small x, 

— V = 2a = (2mc/x)^''^ , m = , (46) 

with coresponding field components 



(47) 



In this case, although the magnetic field grows with decreasing distance from the origin, the 
magnetic force remains negligible in comparison with the gravitational pull of the central 
point mass. Equations (^) and (^Tj) also demonstrate the interesting fact that the local 
(differential) mass-to- flux ratio A = a/bz oc x^/^ — as x — 0^. However, the ratio of the 
spatially integrated mass and flux diverges as the origin is approached, consistent with the 
physical picture of mass accumulating at the center after decoupling from the magnetic 
field. 

Under the fiux-freezing conditions of ideal MHD, A(x) = const, (see §3.3). It may at 
first seem counterintuitive that, in the nonideal case, the local mass-to-fiux ratio decreases 
with decreasing x, since one could expect ambipolar diffusion to become progressively more 
efficient at loading magnetic field lines with matter as one approached the densest parts of 
the core near the center. Indeed, this is what happens before point-mass formation (e.g., 
CM94; CM95). However, after PMF the center acts as a sink of matter, so the magnetic 
field lines become increasingly depleted of mass as one gets closer to the origin (although 
the integrated mass-to-fiux ratio m/ip does increase with decreasing x). 

Our numerical procedure for obtaining the solution in this case has been to integrate 
the differential equations inward and outward from an intermediate point near the origin, 
connecting them to the free-fall solution (eqs. and |^) at small x and to the pivotal 



state (eq. |3^) at large x. Figures 2e-2h show the derived solution for the parameter set 
Ao = 2.9, —Vo = 1.0c, ^_i7 = 5, and rho = 3. These parameter values are typical of the 
outer layers of a supercritical core in the CK simulations and are consistent with physical 
quantities deduced from observations of protostellar cores (e.g., Crutcher et al. 1993, 1994). 
It is seen that, while the solution is everywhere continuous, it exhibits an abrupt change 
in the flow parameters at a certain value of x. This transition, which is characterized by 
a continuous evolution of the physical parameters, a neutral flow velocity that remains 
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supersonic throughout, and a faster deceleration of the ions than of the neutrals, represents 
a C-shock (e.g., Draine & McKee 1993, Smith & Mac Low 1997). A similar ambipolar 
diffusion-induced C-shock has been found in the numerical simulations of CK (see §4 for a 
detailed comparison between the self-similar solution and the CK results). We identify the 
position of the shock with the location of the abrupt drop in the ion velocity and determine 
the shock speed from the motion of this transition relative to the origin. The numerical 
integration yielded Xgh = Vsh = 0.3 and rhc = 5.9, or, equivalently. 



As an aid in visualizing this self-similar solution, we present in Figure 3 time sequences of 
the disk density and radial velocity during the collapse up to the time when a 1.1 Mq star 
is assembled at the center. 

It is worth emphasizing that the C-shock we have found is required by ambipolar 
diffusion and is not the shock discussed in LS (which, as we pointed out in §3.3, is avoided 
altogether in an isopedic disk that has a nonzero initial infall velocity). In the present 
(diffusive) case the magnetic field cannot remain attached to the matter as mass starts 
to accumulate in the origin, so the field decouples from the gas as soon as a point mass 
appears at the center.]^ The region of decoupled flux grows outward at a fraction of the 
speed of sound and a C-shock develops at its outer boundary. An interesting and somewhat 
counterintuitive result is that the accretion rate onto the central point mass is smaller 
(albeit only slightly) in the case with ambipolar diffusion than in the fiux-frozen case 
(compare Figs. 2c? and 2h). The reason is that, as we have noted, the local mass-to-fiux 
ratio decreases when the magnetic field decouples from the matter. The corresponding 
increase in the magnetic force slows down the inflow and thereby reduces the inflow rate. 
Note that the two components of the magnetic force, the field tension and the field pressure 
gradient, are everywhere of comparable magnitude, so one cannot neglect either of them in 
the analysis. 



^^This behavior can also be understood by comparing the ambipolar diffusion and gravitational contraction 
timescales, whose ratio after PMF decreases with diminishing distance from the origin — signaling the onset 
of dynamically important ambipolar diffusion on small scales — on account of the decrease in the free-fall 
time that is brought about by the mass accumulation at the center (see §1 and §3.3 in CK). 





3/2 



Mq per lO*' yr . 
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4. Discussion 

It is of interest to compare our self-similar model results with the detailed numerical 
simulations of CK. The model they present is of a ~ 5 Mq molecular cloud core 
with temperature T = 10 K. Just prior to PMF {t 0~) the accretion rate is 
M = 5.4 Mq per 10^ yr and the mass-to-flux ratio at the center is A = 3.6 — note, however, 
that unlike the self-similar model, neither M nor A are uniform at that time in the model 
core of CK (see their Figs. Ic and Ih). Immediately after PMF the accretion rate rises to 
9.4 Mq per 10^ yr, and it then decreases to 5.6 Mq per 10^ yr by the time a 1 Mq protostar 
forms at the origin (which occurs at t = 1.5 x 10^ yr). Rapid ambipolar diffusion in the 
inner flux tubes halts the inward advection of magnetic flux, which piles up and propagates 
outward as a hydromagnetic disturbance. As the front of piled-up flux moves out to larger 
radii, the magnetic field behind the disturbance and the neutral-ion coUisional coupling 
become strong enough to affect the infalling neutrals, and a shock forms in the core. The 
shock is of the C type (e.g., MuUan 1971; Draine 1980): the infall speed of the neutrals 
(in the shock frame) remains supersonic while the infall speed of the ions is much smaller 
than the ion Alfven speed. By the time a 1 Mq protostar has formed at the origin in their 
typical model calculation, the shock is located at a radius ~ 3.5 x 10'^ AU and propagates 
outward with a speed (in the stellar rest frame) ~ 0.7 C = 0.13 km s^^. (For comparison, 
eq. 1^ yields Mc ~ 1.4 Mq, ~ 1-8 x 10^ AU, and Ugh ~ 0.06 km s"-^, respectively. 



for the accumulated central mass, shock location, and shock speed at that time.) As the 
hydromagnetic disturbance propagates outward, the postshock accretion rate decreases due 
to the fact that the neutrals are "hung up" (i.e., decelerated) in the region of amplified 
magnetic field behind the shock front (see Fig. 6g of CK). This is also what happens in the 
self-similar model (Fig. 2h). As discussed in CK (see also Li & McKee 1996), the postshock 
region is potentially susceptible to interchange instabilities (i.e., the gravitationally-induced 
interchange of magnetic flux tubes; e.g.. Spruit & Taam 1990; Spruit, Stehle, & Papaloizou 
1995): this issue, however, cannot be fully addressed until nonaxisymmetric collapse 
simulations are performed. 

In the limit r (or x 0) the column density and the infall speed of the neutrals 
are both oc r^/^ (oc x^/^, see eq. and Figs. 2e, 2/, 3a, and 36), reflecting the "undiluted" 
free-fall collapse induced by the central point mass. A similar behavior is also revealed in 
the numerical simulations (see Figs. 6c and 6e of CK). 



^^This time is equal to 1.1 Tgr, where Tgr — [t^.i^/GMq]^^'^ is the gravitational contraction (« free-fall) 
time at the location within the core, which contains 1 Mq at t = 0. This result is similar to that 
observed in the nonmagnetic collapse solutions of Shu (1977), Hunter (1977), and Foster & Chevalier (1993). 
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As mentioned in §3.4, ambipolar diffusion "unloads" the mass in the effectively 
stationary (see Fig. 2e) ffux tubes behind the shock front. As a result, the local mass-to-flux 
ratio A(x) = a{x)/hz{x) is no longer a monotonically decreasing function behind the shock 
(see Figs. 2/, 2g, and 36); this behavior is also seen in the CK simulations (see their Fig. 
8a). 

One obvious advantage of the simple self-similar model is that one can scan the space 
of model parameters much more easily than is possible in the case of the time-consuming 
numerical simulations. For instance, the typical model presented in CK took ~ two weeks 
of computer time (running in background on an SGI R-4000 Indigo workstation) for the 
central mass to accumulate IMq of material. By contrast, calculation of our self-similar 
models is usually completed within half a minute (for typical models) when run on the same 
workstation. 

We note that Li & McKee (1996) suggested that Ohmic dissipation [a process that, 
for the assumed cloud composition and cosmic ray ionization rate, becomes important for 
densities n ^ 10^^ cm~^ (e.g., Nakano & Umebayashi 1986a,b), or, equivalently, on scales 
r ^ 1 AU] would halt the accretion of flux onto a protostar, with the flux consequently 
presenting an obstacle to the neutrals and thereby giving rise to a shock. However, as 
shown by CK and verified by our self-similar model, ambipolar diffusion is sufficiently 
efficient at much lower densities (or much larger radii) to halt flux advection and cause 
the formation of a hydromagnetic shock, independent of the effect of Ohmic dissipation. 
Despite the misidentification of the shock formation mechanism and the fact that the 
induction equation was not solved for the magnetic field structure, the simplified model of 
Li & McKee was found by CK to provide a reasonably good approximation to the results 
of their numerical simulations. 

As we have already noted in §3.1, the numerical results of core collapse and point-mass 
formation presented in CK can be described by a self-similar model only in an approximate 
sense. In particular, the condition of the spatial uniformity of u, M, and A at the instant 
of PMF, as well as the scaling nj oc n^/^, are generally not strictly valid in the cores 
of the models presented in CM94, CM95, BM94-BM95b, and CK. Furthermore, the 
temporal behavior of the physical quantities in our self-similar solution deviates from 
that in CK. For example, CK find that in the region of free-fall collapse near the central 
point mass —0.76 < din a /dint < —0.55, which is different from the self-similar result 
din a /dint = —1/2 (see eqs. and [0). It is worth stressing, however, that, in spite 
of these differences, the self-similar solution accurately describes the formation of an 
outward-propagating ambipolar-diffusion shock as well as the other basic characteristics of 
the evolution of a protostellar core during the post-PMF epoch. As discussed in CK, the 
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numerical simulations of this model are consistent with observations of accreting protostars. 
In particular, the inferred evolution of the mass accretion rate is similar to that deduced 
for the protostellar objects HL Tauri and L1551-IRS5, and the calculated magnetic field 
structure agrees with polarimetric observations of dense protostellar cores (see § 4.2 of CK). 

5. Conclusion 

In this paper we have presented a self-similar solution of the collapse of a magnetized 
molccTilar cloud core (assumed to also be nonrotating and isothermal) that, for the first 
time, incorporated the effects of ambipolar diffusion in a self-consistent manner. We have 
focused on the post-PMF (point-mass formation) phase of the collapse of a disk-like core, 
noting that Basu (1997) had previously explored the self-similar nature of the collapse 
before a central mass (i.e., a protostar) first appears at the origin. We clarified the 
distinction between the ideal and nonideal MHD cases by plotting the singular lines in 
the position-velocity space and showing that they correspond to different critical speeds 
(the magnetosonic speed and thermal sound speed in the ideal and nonideal problems, 
respectively). We obtained a solution for the ideal (flux- frozen) case that exhibits a 
split-monopole field topology near the center. This solution differs from the one obtained by 
Li & Shu (1997) in that it involves no shocks. We showed that the shock in the LS solution 
is a direct consequence of their assumption that the core at the time of PMF is described 
by a stationary density distribution (corresponding to a singular isothermal toroid), and we 
pointed out that a shock will generally not be present under the more realistic assumption 
of a nonzero inflow speed near the origin at that instant. We demonstrated, however, that 
a shock is a generic feature of the solution in the nonideal (ambipolar diffusion) case. This 
(C-type) shock is a direct consequence of the action of ambipolar diffusion in the central 
region of the core following PMF: the magnetic diffusivity decouples the fleld from the 
matter, causing the gas to free-fall to the center (where it accumulates in a point mass) 
and the field to stay behind and drive a shock outward. We have compared this solution 
with the results of the numerical simulations of Ciolek & Konigl (1998) and confirmed that, 
while the more realistic numerical models are not strictly self-similar, our simplified solution 
nevertheless captures the main features of the core evolution after PMF. 

This work was supported in part by NASA grants NAG 5-2766 and NAG 5-3687. 
Helpful comments by the referee. Prudence Foster, are gratefully acknowledged. 
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Figure Captions 

Fig. 1. — A plot of the function TZ{X), which determines how different disk radii are 
weighted in the expressions for the gravitational acceleration (eq. [^) and the surface 
radial magnetic field (eq. |Il8|). The function is normalized by /g°°7^(X)dX = 1. 

Fig. 2. — Self-similar collapse solutions for an ideal-MHD (r/ = 0) disk (a) — {d) 
and for an ambipolar diffusion-dominated (r/ = 0.3) disk (e) — (/i), for initial conditions 
represented by the parameters Vo = —1.0, Aq = 2.9, and rho = 3. Plotted as a function 
of the similarity variable x are the radial speeds (in units of the speed of sound C) of the 
neutrals {v, solid curve) and ions (fj, dashed curve) in panels (a) and (e), the dimensionless 
surface density a in panels (b) and (/), the normalized vertical {b^, solid curve) and radial 
surface (&r,s, dashed curve) magnetic field components in panels (c) and (g), and the mass 
accretion rate into the center (M, for a temperature T = 10 K) in panels {d) and {h) . Note 
that the ions comove with the neutrals in the ideal-MHD case (panel [a]). The singular 
curves {open circles) in panels (a) and (e) correspond to straight lines in nonlogarithmic 
units [x — V = {1 + 2A~^)~^ in the ideal-MHD case and a; — f = 1 in the presence of 
ambipolar diffusion]. Both solutions are continuous and have supersonic neutral speeds for 
all X. The abrupt ion deceleration in the nonideal solution gives rise to a strong C-shock. 
No shock appears in the ideal solution. At large x (equivalently, large r), v —>■ const, 
and a oc 6^ ~ br^g oc x~^. At small x (equivalently, small r), the diffusive solution differs 
from the ideal-MHD one. Both solutions contain a point mass at the center and therefore 
both exhibit free-fall-type profiles (f oc a oc near the origin. However, the nonideal 

solution has bz ~ &r,s oc x~^, whereas the ideal-MHD one {bz oc a) contains a split magnetic 
monopole at the center (bz oc and fe^.s oc 

Fig. 3. — The distributions of the core infall speed u and surface number density 
at times t = 0, 2, 4, 6, 8, 10, and 12 x 10^ yr after a point mass forms at the center. By the 
end of the displayed evolution, 1.1 solar masses have accumulated at the origin. 
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